Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 



Printed 1 February 2008 



(MN M^K style file v2.2) 



Non-linear radial oscillations of neutron stars: 
Mode- coupling results 



U. Sperhake^, P. Papadopoulos^, N. Andersson^ 

(1) Department of Mathematics, University of Southampton, Southampton S017 IB J, United Kingdom 

(2) School of Computer Science and Maths, University of Portsmouth, Portsmouth POl 2EG, United Kingdom 



Accepted. Received; in original form October 2001 



ABSTRACT 

The non-linear behaviour of oscillation modes in compact stars is a topic 
of considerable current interest. Accurate numerical studies of such phe- 
nomena are likely to require powerful new approaches to both fluid and 
spacctime computations. We propose that a key ingredient of such meth- 
ods will be the non-linear evolution of deviations from the background 
stationary equilibrium star. We investigate the feasibility of this approach 
by applying it to non-linear radial oscillations of a neutron star, and explore 
numerically various non-linear features of this problem, for a large range 
of amplitudes. Quadratic and higher order mode coupling and non-linear 
transfer of energy is demonstrated and analysed in detail. 



1 INTRODUCTION 



of attention in the past. Although most results have 



Since the first realisation that the Cepheids are stars 
that pulsate radially, stellar oscillations have provided 
a fertile ground for astrophysicists. Observations ob- 
tained with much improved sensitivity have provided 
a wealth of relevant data in recent years. This pro- 
vides an important theoretical challenge given that the 



nature of the various eigenmodes ( Chanmugan 1977 


Glass & Lindblom 1983J; 


Vath & Chanmugam 1992 


Kokkotas & Ruoff 2001 


), there have also been at- 



tempts to study non-linear features in a perturbative 
way (including quadratic and cubic coupling t e rms) 
flDziembowski 1982t |Perdang fc Blachcr 1983); |Pcr- 



gathered information can be matched against detailed 
stellar models. With the likely advent ot gravitational- 
wave astronomy in the next few years, relativists are 
considering whether a similar program may be feasible 
for compact stars. This is an exciting idea since the 
detection of gravitational waves from a pulsating star 
may shed light on the nature of the equation of state 
at supra-nuclear densities. Although stellar oscillation 
theory has been an active field of research for many 
decades (in particular in the context of Newtonian 
gravity) and there are several monographs covering 
the main results, many crucial questions remain open. 
The main uncertainties concern the behaviour in the 
non-linear regime, e.g., the coupling between different 
modes, the formation of shocks etcetera. The purpose 
of the present paper is to demonstrate the accuracy of 
a new approach to the study of non-linear stellar oscil- 
lations. We apply the new method, in which the main 
focus is on non-linear deviations from the background 
stationary equilibrium star, to radial oscillations of 
neutron stars. This is a problem which, given its sig- 
nificance for the stability of the star, has received a lot 



dang 1983|; |Wcntzcl 19871 [Kumar fc Goldreich 1989 



van Hoolst 1996). [It is also relevant to mention the 
recent application of thi s approach to no n-linear ef- 
fects on inertial modes (shenk et al 2001) as well as 



the fully nonli near s imulations by Font, Stergioulas 
and Kokkotas (2000).] Our new approach provides a 



powerful complement to these studies. We investigate 
the main features that appear in the weak to mildly 
non-linear regimes. This leads to some interesting new 
results regarding non-linear mode-coupling and sheds 
light on two of the main questions in this area: i) What 
is the amplitude at which large amplitude modes sat- 
urate?, and ii) How reliable are expansion methods 
beyond quadratic order in the amplitude? We believe 
our paper provides the first detailed study of these 
effects in full general relativity. 



2 NON-LINEAR PERTURBATIONS 

We model the neutron star as a single component per- 
fect fluid at zero temperature which obeys a polytropic 
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equations of state P = Kp 1 , where K and 7 are con- 
stants. For such a fluid the energy momentum tensor 
is given by T M " = (p + P)u^u + Pg<* v , where u M is 



the four velocity, normalised as u^u^ 



-1. We re- 



strict our consideration to spherical stars undergoing 
radial motions, in which case the four velocity is given 
by w M = [v(t, r), w(t, r), 0, 0]. Furthermore, in radial 
gauge and polar slicing the spherically symmetric line 
element is, 



ds 2 = -X 2 dt 2 + pdr 2 + r 2 (d9 2 + sin C J J 



(1) 

i.e. it depends only on two functions X(t, r) and p(t, r). 
With those assumptions the Einstein field equations 
G^v = 87rT MJ/ and the equations of hydrodynamics 
V M T M ^ = lead to four independent equations 



-i + 4nrp 2 [P+{p + P) p 2 w 2 



^ + Anrp 2 [p + (p + P)p 2 w 2 



p,t + aiip.r + ai2i",r = 61, 



+ Qt2\p.r + OllW.r = 62 ■ 



(2) 



(3) 
(4) 



(5) 



Here an, «i2, 021, 61 and 62 (see Sperhake (2001) 



for their exact form) are functions of the fundamen- 
tal variables A, p, p and w, but not of their deriva- 
tives, so that Eqs. (||), (|5|) form a quasi- linear system. 
The solution of Eqs. (0-(^|) requires suitable bound- 
ary conditions. At the centre of the star we require 
that p, = 1 in order to avoid a conical singularity. 
The velocity w vanishes at the origin due to spheri- 
cal symmetry. At the stellar surface we fix the lapse 
function by matching the line element to an exterior 
Schwarzschild metric which implies that Xp = 1. We 
use typical parameters for the neutron star model: 
K — 150 km 2 and 7 = 2. We set the central density 
of the equilibrium configuration of the neutron star 
to p c = 0.001224 km -2 , which leads to a compactness 
(mass/radius ratio) M/R — 0.19. 

In our analysis of non-linear mode coupling we ap- 
proach radial oscillations of neutron stars satisfying 
equations (^)-(^) in the following manner. We decom- 
pose the time dependent quantities A, p, p and P into 
static background contributions and time dependent 
perturbations according to 



f(t,r) = f(r) + Sf(t,r). 



(6) 



The background quantities are determined by the 
Tolman-Oppenheimer-Volkoff equations, the static 
analogues of Eqs. (^)-(^). We use the background 
equations to eliminate all zero order terms from the re- 
sulting perturbative dynamic equations. The resultant 
equations are equivalent to the original system (^)- 
(^). In particular, the two evolution equations form a 
quasilinear system for 8p and w. By eliminating terms 
of order zero in the perturbations we obtain numerical 



accuracy that is determined by the amplitude of the 
perturbation rather than the static background. This 
is the key advantage of this new approach. All prelim- 
inary tests have verified that this non-linear pertur- 
bation scheme provides much enhanced accuracy over 
a large range of amplitudes. Full details of the new 
method as well as the numerical code and its calibra- 



tion are provided by Sperhake (2001) 



Normally the surface of the star is defined by the van- 
ishing of the pressure, which in the polytropic case is 
equivalent to p = 0. However, if one is using an Eule- 
rian framework and the surface of the star is allowed 
to move, the outer grid boundary does not coincide 
with the surface of the star and this condition cannot 
be applie d easil y. As has been discussed in detail by 
Sperhake (2001) this leads to severe numerical difficul- 
ties, and could trigger artificial shock formation in the 
surface region. In the present study we want to focus 
on the non-linear coupling between various oscillation 
modes. In order to isolate this effect (and avoid any 
artificial effects due to the surface of the star) we use 
a fixed (rather than free) boundary condition, i.e. we 
require w = at the surface. Furthermore we do not 
evolve the low density layers of the neutron star, in 
order to avoid negative total energy densities. The re- 
sulting neutron star model contains about 90 % of the 
mass of the original model. 

The eigenmodes of a dynamic spherically symmetric 
neutron star are described by the linearised version 
of the dynamic equations (^)-(^). It is a well known 
result that the linearised equations lead to a self ad- 
joined eigenvalue problem in terms of the rescaled dis- 
placement vector (,, which is related to our variables 
cf. chapter 26 in Misner, Thorne and 



by C,t = ' 

Wheeler (1973). The solutions of the eigenvalue prob- 



lem form a complete orthonormal system and hence 
we can expand the time dependent £(t, r) resulting 
from fully non-linear evolution in a series of the linear 
eigenmodes 



C(i,r) = ^4,(i)(,(r), 



(7) 



where the coefficients Ai are given in terms of the 
inner product 



Mt)= f R (P + P )^0at,r)Ur)dr. 
Jo r 



(8) 



In Fig.[l] we show the four lowest eigenmodes in the 
velocity field for our truncated stellar model. The ex- 
pansion (Q) provides a useful diagnostic which allows 
us to assess the level of non-linear coupling through- 
out a numerical evolution. We measure the presence 
of mode i at any given time during the evolution by 
calculating the coefficient Ai(t) via numerical integra- 
tion. 
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Figure 1. The profiles of the lowest four eigenmodes in the 
velocity w for the neutron star model with fixed surface. 

3 RESULTS 

In order to study non-linear mode-coupling we evolve 
the nonlinear perturbation equations from initial data 
corresponding to the velocity field Wj of a single lin- 
ear eigenmode (of order j) with amplitude Kj = 
max|A£j(r)/r 2 |, which represents the maximum radial 
displacement of a fluid element inside the star. The 
initial density perturbation Sp is set to zero, while the 
initial values for the metric variables follow from the 
constraints 

3.1 Exciting the fundamental mode 

We first consider the case when the initial data cor- 
respond to the fundamental radial eigenmode. Hav- 
ing evolved this data, we measure the maximal co- 
efficients Ai = max\Ai{£)\ obtained over an integra- 
tion time corresponding to many times the dynamical 
timescale. In Fig.^ we show the maximal coefficients 
obtained for the lowest 10 eigenmodes for excitation 
amplitudes ranging between 1 cm and 50 m. We ob- 
serve weak mode-coupling throughout most of this do- 
main. The fundamental mode itself (Ai), is seen to 
grow more or less linearly with the initial amplitude 
K\, which indicates the absence of significant self in- 
teraction. Meanwhile, for higher order modes we can 
clearly identify two different regimes: For amplitudes 
below 10 m, all coefficients A2, . . . , A10 grow quadrat- 
ically with the excitation amplitude K\ . At larger am- 
plitudes all eigenmode coefficients except for A2 show 
a transition to power laws with larger index. We have 
illustrated this behaviour in Fig.|2j by modelling the 
coefficients A2, A3 and A4 as power series expansions 
in K\ according to 

Ai = aK? + diKt, (9) 

(the Einstein summation convention is not use here 
or in similar expressions below) where d = {3.6 ■ 
10~ 7 , 3.4-10" 8 , 1.0-10- 8 }anddi = {0, 9.7-1Q- 10 , 1.2- 



A, 

1(T Z 
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Figure 2. The coupling of the fundamental mode to higher 
order modes is illustrated by plotting the coefficients Ai as 
functions of the initial amplitude K\ . For low amplitudes 
the results are well fitted by quadratic power laws while 
higher powers are required in the mildly non-linear regime 
(as indicated by the solid lines). 

10 -11 } for i — 2,3,4. The higher order power laws 
have been obtained by least square fits to the coef- 
ficients Ai after subtracting the quadratic contribu- 
tions CiKi. For the modes i = 5 — 10 the contribution 
of the higher order power law is rather weak which 
makes it difficult to obtain accurate measurements of 
the corresponding exponents. The steepening of the 
curves is, however, still obvious in the figure. It is also 
clear, since the curves for the eigenmode coefficients 
Ai do not intersect in Fig. ^, that the coupling strength 
decreases with the order of the eigenmodes over the 
whole range of amplitudes. 

3.2 Exciting higher order modes 

Next we set initial data in the form of second and 
higher eigenmodes. In this case we still observe the two 
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Figure 3. The eigenmode coefficients obtained for initial 
data in the form of the second eigenmode. We have fit- 
ted the data in the intermediate amplitude range between 
K\ = 1 m and 10 m by quadratic power laws. The top 



-i/2 



panel shows the resulting fits to Cj • K 2 + di ■ K. 



regimes mentioned above. In the weakly non-linear 
regime the eigenmode coefficients are well approxi- 
mated by quadratic power laws. The main difference 
from the previous case is that the second mode (of 
order i) shows a preference to couple to modes of or- 
der 2i. This is particularly clear in Fig. || where the 
coefficients Ae, As and Aio show a much steeper in- 
crease for large amplitudes. For these modes we also 
obtain excellent fits with combinations of power laws 
according to 



A t = aKi + d t K* /2 



2.5-10" 

-9 



4.9-10" 

-it 



(10) 
1.8- 10" 9 } 



where a = {2.9-10" 

and di = {0, 8.7 • 10" 9 , 6.2 ■ 10" 11 , 4.9 ■ 10""} for 
i — 4,6, 8, 10. These fits are shown in Fig.||. 
These results are generally confirmed if the initial data 
is given in the form of the third velocity mode. The 



only difference is that the preferred modes in the mod- 
erately non- linear regime are now those of order 3i. 



4 DISCUSSION 

In this paper we have applied a new non-linear ap- 
proach to the study of stellar pulsation, and studied 
mode-coupling due to non-linear effects by evolving 
initial data corresponding to a single linear eigenmode 
with varying amplitude. Concerning the transfer of 
energy to other modes we have found two distinct 
regimes, a weakly non-linear regime where the excita- 
tion of modes grows quadratically with the initial am- 
plitude and a moderately non-linear regime, which can 
be reasonably well described by power laws of higher 
order. 

The results for the weakly non-linear regime agree 
qualitatively with Newtonian perturbative studies. In 
the analytic study of non-linear mode coupling one 
normally views the eigenmode coefficients as har- 
monic oscillators and the non-linear interaction be- 
tween eigenmodes is represented in the form of driv- 
ing terms which are quadratic or of higher order in 
the amplitudes (see for example van Hoolst (1996)) 



d 2 A t 
dt 2 



+ ujfAi = cfAjAk + d J iAjA k Ai + 



(11) 



where the <£ , d 



■jkl 



are the quadratic, cubic and 



higher order coupling coefficients and summation over 
j, k, I is assumed. In our analysis the initial data con- 
sists of one isolated eigenmode j, so that the right 
hand side can be approximated by c%A 2 + diAj + . . . 
In analytic studies this series expansion is normally 
truncated at second or third order. The omission of 
higher order terms is justified in the weakly non-linear 
regime, where our fully non-linear simulations confirm 
that quadratic terms in the initial amplitude dominate 
the coupling between eigenmodes. This is no longer 
true, however, in the moderately non-linear regime, 
where higher order terms are more important. Our re- 
sults allow us to define the transition to this regime, 
manifested by the breaks in the curves in Figures H 
and |j| As is clear from the figures, the moderately 
non-linear regime corresponds to initial mode ampli- 
tudes above 10 m or so. We note that the correspond- 
ing Mach number is of the order of 0.01. This agrees 
well with investigations of Newtonian stars (Kumar 
& Goldreich 1989) which assume a Mach number of 
0.1 as the limit of applicability of semi-analytic mode 
coupling methods. 

Furthermore, our results indicate that the non-linear 
couplings would be poorly captured by polynomial ex- 
pansions in the mildly non-linear regime. We observe 
significant excitation of higher order modes (eighth, 
tenth etc) which can only emerge from very high-order 
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couplings. To quantify the associated coupling coef- 
ficients in a perturbative calculation would be very 
difficult. 

We have also observed that, given an initial mode j, 
the coupling to modes nj is particularly efficient in the 
moderately non-linear regime. This is naturally inter- 
preted as a resonance effect. In analogy with the sim- 
ple problem of a single forced oscillator we can assume 
that resonance occurs for any mode whose frequency 
is an integer multiple of the driving frequency ujj in 
the general non-linear case, i.e. we can schematically 
write the eigenmode coefficients in the form 

Mt)~y\ 2 F ? » . (is) 

where the F n will depend o n th e frequencies (cf. 
Eqs. (18), (19) of van Hoolst fll996| ). In our case the 
external force is provided by the non-linear coupling 
to the initial mode j, as indicated in ([l2|). We therefore 
obtain resonance if u>i — nujj . For our simple neutron 
star model, the eigenfrequencies of radial neutron star 
oscillations are almost equally spaced and we can use 
ujj w (jtOi)/i for i,j > 2 as a reasonable approxima- 
tion. The condition for resonance then becomes i = nj 
which is exactly what we have observed. 
As one of the main results of this paper, we emphasise 
the new perturbative approach that enabled us to ob- 
tain highly accurate, fully non-linear, evolutions over 
a large range of amplitudes. This a pproach was dis- 
cussed in detail by Sperhake (2001), and we plan to 
present further results regarding mode-coupling and 
non-linear shock formation in future papers. In prin- 
ciple, this technique can be applied to any physical 
problem that involves a non-trivial stationary limit 
and we expect it to prove a valuable tool in many non- 
linear problems. We note that the accuracy improve- 
ments are independent of the numerical discretization 
used (here, second order, centred finite differences). In 
combination with methods suitable for smooth oscilla- 



tory solutions (Gourgoulhon 1991; Gourgoulhon et al 
1995| ), we would expect a dramatic expansion of the 
applicability of non-linear simulations to relativistic 
stellar pulsations. 

A problem for which our new approach may prove use- 
ful concerns unstable modes of rotating neu tron s tars 
(eg. the r-modes, see Andersson & Kokkotas (2001) for 



a review) . One of the most important questions raised 
in connection with the r-modes concerns the ampli- 
tude at which an unstable mode saturates. Recently, 
direct three dimensional numerical simulations have 



been brou g ht to bear on this prob lem (Stergioulas 



Font 2001 



Lindblom et al 2001). The picture that 



emerges from these studies is, however, not conclusive. 
Both studies suggest that an unstable r-mode satu- 
rates at an extremely large amplitude (corresponding 
to waves of a height of several hundred meters in a 
star spinning near the breakup limit) due to shocks 



forming in the surface region. That such "wave break- 
ing" would occur once a mode reaches a large ampli- 
tude is likely, but simulations must isolate the true 
physical behaviour near the surface of a star, from 
numerical artifacts associated with the rapid decrease 
i n the density (for a detailed discussion see Sperhake 
(2001)). In fact, it is interesting to contrast these re- 



sults with those of the present work that suggest that 
nonlinear effects are highly relevant already at wave 
amplitudes of order 10 m. We believe that a suitable 
generalisation of the method used in this paper, that 
provides unprecedented accuracy for a large range of 
wave amplitudes, could prove extremely useful for the 
study of unstable non-axisymmetric modes and plan 
to address such problems in the near future. 
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